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■i  i  1?  ae’.xlopniont  of  nurnerica!  methods  for  solving  problems  in  inateriids  science  is  a  subject  which, 
v.  itii  ii  ;  /  exceptions,  is  in  its  infancy.  This  work  covers  one  of  the  basic  computational  problems  in 
-  ' e  hold  which  is  that  of  devising  algorithms  for  computing  microstructures.  Our  efforts  have  all 
boeii  diii'cled  towards  computing  the  particular  microstructures  associated  with  martensitic  phase 
iransitionc.  Within  this  framework  have  encountered  both  algorithmic  and  analytical  problems, 
a  number  of  which  have  been  satisfactorily  solved. 


Our  work  addresses  static  problems.  W'e  are  not  currently  able  to  handle  fuU  three  dimensional  time 
dependent  cases.  The  typical  stationary  problem  calls  for  the  minimization  of  a  highly  nonconvex 
bulk  energy  functional.  The  extremizing  functions  are  vector  valued,  representing  displacement 
fields.  The  construction  of  the  bulk  energy  functionals  is  a  specialized  task,  depending  on  a  detailed 
knowledge  of  the  properties  of  materials.  The  functionals  used  in  our  work  were  derived  by  R.  James, 
based  on  ideas  of  J.  Ericksen.  The  functionals  contain  information  about  the  structure  of  the  crystal 
lattice  and  several  material  specific  parameters.  These  functionals  are  relatively  complex,  being, 
for  instance,  eighth  order  polynomials  in  the  displacements. 

The  research  supported  in  the  grant  has  addressed  two  main  issues.  The  first  issue  concerns 
techniques  for  discretization  of  the  bulk  energy  functionals  and  the  second  concerns  algorithmic 
issues  relevant  to  nonconvex  optimization.  In  this  we  discovered  some  unusual  and  apparently 
paradoxical  results.  We  will  describe  these  results  now. 

As  reported  in  [1]  the  initial  discretization  used  biUnear  finite  elements  in  a  planar  model.  In  addi¬ 
tion  to  the  finite  element  approximation  it  is  necessary  to  use  some  form  of  numerical  integration 
to  evaluate  the  energy  since,  as  mentioned  above  it  is  a  polynomial  of  high  degree.  Following  the 
work  of  others,  we  used  a  midpoint  quadrature  approach  combined  with  bilinear  elements  on  a 
uniform  square  mesh.  Minimizing  the  resulting  function  of  the  mesh  variables  (several  thousand  of 
them)  leads  in  principle  to  the  desired  microstructure.  Our  important  discovery  here  was  that  the 
‘microstructures’  computed  by  others  using  this  scheme  are  not  microstructures  at  all,  but  are  rep¬ 
resentations  of  numerical  instabilities  of  the  discretization  technique.  This  situation  occurs  because 
both  the  microstructure  and  the  numerical  instability  are  manifested  as  grid  scale  oscillations.  It 
can  be  hard  to  distinguish  the  two  in  many  cases,  but  we  constructed  an  example  where  there  can 
be  no  doubt  that  the  computed  result  represents  the  eifect  of  numerical  instability,  although  it  has 
the  physical  appearance  of  a  typical  microstructure.  The  report  [1]  contains  an  explanation  of  the 
source  of  this  instability  and  shows  how  to  avoid  it  it  certain  cases.  It  seems  unlikely  that  it  can 
be  avoided  in  general.  If  this  is  true  it  means  that  the  original  scheme  is  seriously  flawed  and  v/e 
do,  in  fact,  consider  this  to  be  the  case. 

'£o  avoid  the  difficulties  of  the  approach  of  the  previous  paragraph  we  decided  to  eliminate  the 
Integration  errors  entirely  by  the  use  of  linear  finite  elements  on  triangles.  This  introduces  another 
element  of  uncertainty  in  that  -with  triangular  meshes  the  grid  orientation  may  play  a  role  in  the 
computed  microstructure.  Our  experiences  with  this  are  reported  in  [2].  The  major  difficulty  is  that 
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wkm  tlie  mesh  disgonsls  are  not  aligned  with  the  microstructure  they  nevertheless  ue  somehow 
elde  Ip  *cs^Uire*  the  microstmeture  and  produce  a  result  in  which  it  is  aligned  with  the  mesh 
;d||0tNU4s.  This  can  occur  even  when  the  physical  microstmeture  is  meant  to  lie  at  90  degrees  to 
tike  madi  ^fiagonak.  The  smnewhat  subtle  reasons  for  the  occurence  of  this  phenomenon,  explained 
hi  we  to  he  found  in  the  existmice  of  near  minimizers  in  the  energy  sense  which  are  far  from  the 
phydi^  minimiiers.  At  present,  the  only  available  solution  to  this  difficulty  is  to  compute  on  very 
fora  meshes.  Unfortunately,  this  is  an  expensive  remedy.  We  have  developed  a  theory  to  explain 
gWAtitativ^  how  fine  the  mesh  needs  to  be  in  particular  cases  and  the  techniques  for  this  are 
presented  In  [2]. 

hi  the  course  of  p«rfonning  the  computations  of  the  two  previous  paragraphs  a  suprising  effect  of 
numwied  optboization  algorithms  was  noted.  Specifically,  when  applied  to  nonconvex  functionals, 
the  usual  deseeik  algorithms  can  terminate  at  a  saddle  point  at  the  functional  rather  than  at  a  local 
ndnimum  as  would  be  expected.  No  warning  of  this  situation  is  provided.  It  simply  appears  to  the 
nser  that  a  minimum  (local)  has  been  found  whoeas  in  fact  the  computed  solution  is  only  a  saddle 
pOiffit.  We  p>n8tracted  a  simple  model  to  demonstrate  this  effect  and  provided  some  analysis  of 
foe  draraeter  the  saddle  point.  This  work  is  bdng  written  up  for  publication  [3].  We  are  unsure 
how  pervasive  this  phenomemm  is  likdy  to  be  in  general,  but  we  do  have  full  scale  microstructure 
cemptttafoms  in  which  it  occurred. 

difficulties  of  cmnputing  mkrostructure  in  the  form  of  grid  scale  osdllations  are  sufficiently 
that  it  beoesnes  of  interest  to  sedc  ahemative  formulations.  In  one  such  formulation  we 
sefo  averagss  cl  the  osciilations  rather  than  the  osctllations  themselves.  A  technique  for  this  was 
ptpposed  in  our  work  [4}.  We  formulated  an  approach  which  is  suitable  for  directly  computing 
foe  Yoog  measures  whifo  describe  the  mesh  scale  oscillations.  The  calculations  in  this  case  take 
|fokoe  on  a  fixed  ipefo  which  h  mA  unduly  fine  since  the  objects  bedng  computed  are  relatively 
iios%  varying.  O^e  ibawback  this  iq»pioach  is  that  rince  only  averages  are  being  found  the 
ndcroatructure  must  somehow  be  obtained  from  the  averages.  For  design  purposes  however,  the 
'muspgP-  way  be  suffident  information.  We  had  some  success  in  computing  with  this  approach  in 
ffiffibipri  dtuations.  However,  there  are  more  variables  to  keep  track  of  and  the  optimization 
pKoldan  is  affifficult  one. 

^  Ibyelilld  tfok  pmtiy  cwiiputatkmal  issues,  there  is  the  question  of  assessing  the  theoretical  accuracy  of 
tifo  saitttimtB.  The  mkrostructure  problem  is  far  too  difficult  for  existing  techniques, 

fo  fofd/fora  uik^eil^^  mathraatical  theory  is  not  wdl  devfo>ped  for  the  particular  problems  of 
.smdlw-ti  k  unvdse  to  attempt  to  prove  error  estimates  for  these  problems  at  present, 
fotiiad,  wn^tiidded  to  see  ufoat  can  be  said  about  simpkr  modds.  In  one  case,  reported  in  [5]  we 
uMra  lb  otAaln  ifomous  reralts  for  a  omvexified  problem.  Our  estimates  are  norm  estimates 
bad  dak  bettrar  ihaa  est^  hi  weak  topolof^. 

fo  of  tUuciwtisnikm 

1^^  to  dkbrdizatkm  uses  pairs  of  dual  meshes  to  closely  enforce  physical  laws 
tfdu^ue  is  designed  specffically  around  discretizing  .the  vector  fidd 
^  varkms  combinations.  One  lug  advantage  of  covdnme 

fo jfold  ltii^^|d^t}Uas.stt€h  as  div(curi)sO  and  cUrl(grad)aO  hi^  natural  discrete 
iabiti  fonfoi  Jk  sofid  Iheonffical  fnmewofo  has  been  built  for  these  methods  and  they 

mate  a  numbm  of  omtribufoms  to  this  effort.  In  the 
contributioiis  which  are  described  bdow. 

fou  covdume  sfoeme  for  three  dimendonal  div-curl  problems 


litv.  mcah  <Hago35al^  are  aol  aligned  with  the  microstructure  they  nevertheless  are  soriiehow 
to  capture  the  microstructure  and  produce  a  result  in  which  it  is  aligned  with  the  mesh 
u  Cj^oiials.  iliis  can  occur  even  when  the  physical  nucrostructure  is  meant  to  lie  at  90  degrees  to 
i!’.;'  intsii  diugonals.  The  somewhat  subtle  reasons  for  the  occurence  of  this  phenomenon,  explained 
'  ‘  .  to  be  found  in  the  existence  of  near  minimizers  in  the  energy  sense  which  are  far  from  the 

pi;  \  .::cal  lainimizGis.  At  present,  the  only  available  solution  to  this  difliculty  is  to  compute  on  very 
fiiU!  meshes.  Ursfortunately,  this  is  an  expensive  remedy.  We  have  developed  a  theory  to  explain 
(juantitatively  how  fine  the  mesh  needs  to  be  in  particular  cases  and  the  techniques  for  this  are 
presented  in  [2]. 


in  the  course  of  [)erforming  the  computations  of  the  two  previous  paragraphs  a  suprising  effect  of 
numerical  optimization  algorithms  was  noted.  Specifically,  when  applied  to  nonconvex  functionals, 
the  usual  descent  algorithms  can  terminate  at  a  saddle  point  of  the  functional  rather  than  at  a  local 
minimum  as  would  be  expected.  No  warning  of  this  situation  is  provided.  It  simply  appears  to  the 
user  that  a  minimum  (local)  has  been  found  whereas  in  fact  the  computed  solution  is  only  a  saddle 
point.  We  constructed  a  simple  model  to  demonstrate  this  effect  and  provided  some  analysis  of 
the  character  of  the  saddle  point.  This  work  is  being  written  up  for  publication  [3].  We  are  unsure 
how  pervasive  this  phenomenon  is  likely  to  be  in  general,  but  we  do  have  full  scale  microstructure 
computations  in  which  it  occurred. 


The  difficulties  of  computing  microstructure  in  the  form  of  grid  scale  oscillations  are  sufficiently 
great  that  it  becomes  of  interest  to  seek  alternative  formulations.  In  one  such  formulation  we 
seek  averages  of  the  oscillations  rather  than  the  oscillations  themselves.  A  technique  for  this  was 
proposed  in  our  work  [4].  We  fonnnlated  an  approach  which  is  suitable  for  directly  computing 
the  Young  measures  which  describe  the  mesh  scale  oscillations.  The  calculations  in  this  case  take 
place  on  a  fixed  mesh  which  is  not  unduly  fine  since  the  objects  being  computed  arc  relatively 
slowly  varying.  One  drawback  of  this  approach  is  that  since  only  averages  are  being  found  the 
microstructure  must  somehow  be  obtained  from  the  averages.  For  design  purposes  however,  the 
averages  may  be  sufficient  information.  We  had  some  success  in  computing  with  this  approach  in 
several  different  situations.  However,  there  are  more  variables  to  keep  track  of  and  the  optimization 
problem  is  a  difficult  one. 

Beyond  the  purely  computational  issues,  there  is  the  question  of  assessing  the  theoretical  accuracy  of 
the  computed  solutions.  The  full  microstructure  problem  is  far  too  difficult  for  existing  techniques. 
In  fact,  the  underlying  mathematical  theory  is  not  well  developed  for  the  particular  problems  of 
interest  and  so  it  is  unwise  to  attempt  to  prove  error  estimates  for  these  problems  at  present. 
Instead,  we  decided  to  see  what  can  be  smd  about  simpler  models.  In  one  case,  reported  in  [5]  we 
were  able  to  obtain  rigorous  results  for  a  couvexified  problem.  Our  estimates  are  norm  estimates 
and  thus  better  than  estimates  in  weak  topok^es. 


2.  Cofvolume  techniques  of  discretization 

The  covolume  approach  to  discretization  uses  pairs  of  dual  meshes  to  closely  enforce  physical  laws 
at  the  discrete  level.  The  technique  is  designed  specifically  around  discretizing  the  vector  field 
cqterators  div,  grad  and  curl  and  thmr  various  combinations.  One  big  advantage  of  covolume 
discretizations  is  that  vector  identities  such  as  (iiv(curl)=0  and  curl(grad)sO  have  natural  discrete 
analogs.  In  recent  times,  a  soUd  theoretical  framework  has  been  built  for  these  methods  and  they 
have  been  applied  in  several  fields.  We  have  made  a  number  of  contributions  to  this  effort.  In  the 
current  grant  period  we  have  made  the  contributions  which  are  described  bdow. 

In  [6]  we  present  an  error  estimate  for  the  covolume  scheme  for  three  dimensional  div>curl  problems 
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QH  taftarskednl  mashes  is  three  dimoisions.  The  three  dimeasional  div-carl  problem  is  a  notoriously 
4Mk»lt  <me.  It  is  ov^etermined  since  it  has  four  equations  bat  only  three  dimensions.  There  are 
fenTt  if  nay,  discretiinticms  which  deal  with  the  system  directly.  Usually,  potentials  or  Biot-Savart 
hitniralB  are  used.  Another  indirect  approach  can  be  based  on  least  squares  although  not  without 
On  the  other  hand  our  algorithm  provides  a  straightforward  way  to  obtain  a  solution  to 
mh  qrstems  and  the  error  estimate  shows  that  the  discretization  is  effectively  optimal. 

Oae  proUmn  with  covdume  methods  is  that  unlike  finite  elements  we  do  not  have  schemes  of 
arbitrary  wder  of  accura^  for  general  meshes.  As  part  of  an  effort  to  provide  such  schemes  we 
have  devek^Mid  a  higher  oovolume  scheme  for  div-carl  systems  on  r^ular  meshes.  This  scheme 
is  baSt  <m  the  standard  schmne.  Its  accuracy  is  fourth  order  in  the  mesh  spacing  compared  with 
second  order  for  the  standard  scheme.  We  have  a  r^orons  proof  of  the  order  of  the  scheme.  This 
is  repmrted  m  [7].  The  algmithm  is  equally  ^plkable  in  two  and  three  dimensions.  However,  we 
have  BO  yet  succeeded  in  extending  it  to  triangular  and  tetrahedral  meshes  in  any  concise  way.  If 
wme  posdBde,  the  lesnlting  scheme  would  most  likely  be  third  order  as  opposed  to  first  order 
for  the  standard  scheme. 

Am  aeathHied  frevkmsl^,  covohime  schemes  use  dual  mesh  systems.  The  prototypical  such  mesh 
|Wbr  is  a  Voronoi-Ddannay  system,  l^ese  mcdiee  have  numerous  applicaticms  in  computational 
jgieindiy  and  have  been  used  as  finite  donent  meshes  for  a  number  of  years.  Thdr  use  in  the 
eesudsme  scheme  is  associated  with  the  gecnnetrical  property  that  the  edges  of  the  Voronoi  diagram 
>  are  octhogenal  to  the  faces  of  the  Ddaunay  triangulation  (or  tesselation)  and  reciprocally.  The  usual 

/  'SitBaddft  envisagss  n  convex  d<miatn  to  be  triangulated,  but  of  coarse  this  is  insufficient  for  partial 
dSffinrsntial  eqnatkms  where  restrictions  on  geometries  cannot  be  td»ated.  For  this  reason  we  have 
'  developed  n  new  algorithm  for  making  Vturemd-Ddaunay  mesh  systems  for  arbitrary  dmnains.  Our 
|sm|r  ilgcfithin  »  based  on  condating  a  set  of  nodes  to  a  uniform  mesh.  (The  uniform  mesh  is  used 
as  ai(^  ^  forming  the  Voioitoi-Ddaunay  system.  It  has  nothing  to  do  with  the  discretization 
?  ilSiif4  Ws  have  bem  able  to  show  that  our  a|p>rithm  is  optimal  in  time  and  space.  It  is  highly 
sid^dde  fat  imrtial  dffierentid  ^jfficatkm  haring  a  number  of  local  properties.  For  instance,  we 
hdup  alisad.<^  time  how  maiQr  and  which  dements  have  to  be  changed  when  a  new  pdnt  is  added. 
,  fids  Is  repeated  in  {8]. 

flM^lojdc  {8}  whkdi  we  assembled  and  edited  contains  a  survey  article  for  covolume  algorithms 
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fj^-capturing  effects  In  nucrostructure  oomputation8;(with  N.  Walkington  and 
0  be  sd^tted  to  SIMf  Journal  on  Scientific  and  Statistical  Computing) 
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e^hantes  for  Young’s  variatkuial  problem:(with  N.  Walkington)  (to  appear  in 
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>*ji  tetraht'dra!  meshes  iit  three  dimensions.  The  three  dimensional  div-curl  problem  is  a  notoriously 
di'licult  one.  It  is  overdeterrnined  since  it  has  four  equations  but  only  three  dimensions.  There  are 
few,  if  any,  discretizations  which  deal  with  the  system  directly.  Usually,  potentials  or  Biot-Savart 
integrals  are  used.  Another  indirect  approach  can  be  based  on  least  squares  although  not  without 
diffktdty.  On  the  other  hand  our  algorithm  provides  a  straightforward  way  to  obtain  a  solution  to 
siicli  systems  and  the  error  estimate  shows  that  the  discretization  is  effectively  optimal. 

One  problem  with  covolume  methods  is  that  unlike  finite  elements  we  do  not  have  schemes  of 
arbitrary  order  of  accuracy  for  general  meshes.  As  part  of  an  effort  to  provide  such  schemes  we 
have  developed  a  higher  covolume  scheme  for  div-curl  systems  on  regular  meshes.  This  scheme 
is  built  on  the  standard  scheme.  Its  accuracy  is  fourth  order  in  the  mesh  spacing  compared  with 
second  order  for  the  standard  scheme.  We  have  a  rigorous  proof  of  the  order  of  the  schemi. .  This 
is  reported  in  [7].  The  algorithm  is  equally  applicable  in  two  and  three  dimensions.  However,  we 
have  no  yet  succeeded  in  extending  it  to  triangular  and  tetrahedral  meshes  in  any  concise  way.  If 
this  were  possible,  the  resulting  scheme  would  most  likely  be  third  order  as  opposed  to  first  order 
for  the  standard  scheme. 

As  mentioned  previously,  covolume  schemes  use  dual  mesh  systems.  The  prototypical  such  mesh 
pair  is  a  Voronoi-Delaunay  system.  These  meshes  have  numerous  applications  in  computational 
geometry  and  have  been  used  as  finite  element  meshes  for  a  number  of  years.  Their  use  in  the 
covolume  scheme  is  associated  with  the  geometrical  property  that  the  edges  of  the  Voronoi  diagram 
arc  orthogonal  to  the  faces  of  the  Ddannay  triaagnlation  (or  tesselation)  and  reciprocally.  The  usual 
situation  envisages  a  convex  domain  to  be  triangulated,  but  of  course  *:his  is  insufficient  for  partial 
differential  equations  where  restrictions  on  geometries  cannot  be  tederated.  For  this  reason  we  have 
developed  a  new  algorithm  for  making  Voronm-Ddaunay  mesh  systems  for  arbitrary  domains.  Our 
new  al^rithm  is  based  on  correlating  a  set  of  nodes  to  a  uniform  mesh.  (The  uniform  mesh  is  used 
as  a  tool  for  forming  the  Voronoi-Delaunay  system.  It  has  nothing  to  do  with  the  discretization 
itself.)  We  have  been  able  to  show  that  our  algorithm  is  optimal  in  time  and  space.  It  is  highly 
suitable  for  partial  diflferential  application  having  a  number  of  local  properties.  For  instance,  we 
know  ahead  of  time  how  many  and  which  elements  have  to  be  changed  when  a  new  point  is  added. 
This  work  is  reported  in  [8]. 

In  [10]  we  presented  a  comprehensive  covolume  scheme  for  the  compressible  Navier-Stokes  equations. 
The  book  [9]  which  we  assembled  and  edited  contains  a  survey  article  for  covolume  algorithms 
generally. 
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